version 9.0
cap log close
log using log\bs.log , replace

clear
set seed 365476247

cap prog drop runme
prog def runme

local hypothesis = 0

tempfile main bootsave

use dta\Turnout_District_Level, clear 

/* GENERATING QUANTILES */
gen Margin=SMD_margin_pre
tab PR_district, gen(d_PR_district)
keep if year==1921
*****************

di
reg ch_turnout Margin 
global mainbeta = _b[Margin]
global maint = (_b[Margin] - `hypothesis') / _se[Margin]
predict epshat, resid
gen yhat=ch_turnout-epshat

/* also generate "impose the null hypothesis" yhat and residual */
gen temp_y = ch_turnout - Margin * `hypothesis'
reg temp_y 
predict epshat_imposed, resid
gen yhat_imposed=ch_turnout-epshat_imposed
qui replace yhat_imposed = yhat_imposed + Margin * `hypothesis'

sort PR_district year
qui save `main' , replace

**qui by year: keep if _n == 1
**qui summ
**global numyears = r(N)

cap erase `bootsave'
qui postfile bskeep beta_np t_np t_wild using `bootsave' , replace

global bootreps = 999

forvalues b = 1/$bootreps {

/* first do wild bootstrap */
use `main', replace
qui by PR_district: gen temp = uniform()
qui by PR_district: gen pos = (temp[1] < .5)
/*  these two lines (and commendted t-stat below
are if you don't want to impose the null hypothesis*/
/*
gen wildresid = epshat * (2*pos - 1)
gen wildy = yhat + wildresid
*/
gen wildresid = epshat_imposed * (2*pos - 1)
gen wildy = yhat_imposed + wildresid
qui reg wildy Margin, cluster(PR_district) 
local bst_wild = (_b[Margin] - `hypothesis') / _se[Margin]
*local bst_wild = (_b[post_self] - $mainbeta) / _se[post_self]



/* next do nonparametric bootstrap */
bsample , cluster(PR_district) idcluster(newPR_district)
********replace knr=knr-100*cnr+100*newcnr  /* only relevant if doing fixed effect where knr is panel var */
qui reg ch_turnout Margin, cluster(newPR_district) 
local bsbeta = _b[Margin]
local bst = (_b[Margin] - $mainbeta) / _se[Margin]

/* now do condition on X bootstrap */
/* first randomly sort new clusters. then merge residuals onto main data */
**keep epshat newyear selfemployed
**gen newsort = uniform()
**sort newyear
**qui by newyear: replace newsort = newsort[1]
**sort newsort selfemployed
**qui by newsort: gen first = 1 if _n == 1
**gen newnewyear = sum(first)
**drop newyear newsort
**gen year = 1982 + newnewyear
**sort year self
**rename epshat newepshat
**merge year self using `main'
**gen condy = yhat + newepshat
**qui reg condy selfemployed post post_self , cluster(year)
**local bst_cond = (_b[post_self] - $mainbeta) / _se[post_self]

post bskeep (`bsbeta') (`bst') (`bst_wild')
} /* end of bootstrap reps */
qui postclose bskeep



qui drop _all
qui set obs 1
gen t_np = $maint
gen t_wild = $maint
qui append using `bootsave'

qui gen n = .
foreach stat in t_np t_wild {
qui summ `stat'
local bign = r(N)
sort `stat'
qui replace n = _n
qui summ n if abs(`stat' - $maint) < .000001
local myp = r(mean) / `bign'
global pctile_`stat' = 2 * min(`myp',(1-`myp'))
}

global mainp = norm($maint)
global pctile_main = 2 * min($mainp,(1-$mainp))

local myfmt = "%7.5f"

di
di "Number BS reps = $bootreps, Null hypothesis = `hypothesis'"
display "Main beta" _column(13) "main T"	_column(22) "Main %le" ///
	_column(33) "NP BS %le" _column(44) "fixX %le" _column(54) "wild %le"
di 	%6.3f $mainbeta _column(13) %6.3f $maint  _column(23) ///
	`myfmt' $pctile_main _column(34) `myfmt' $pctile_t_np ///
		_column(55) `myfmt' $pctile_t_wild

end


runme
log close
